Interaction QTL x Year with QTLRel | SS
Interaction QTL x Year with QTLRel | SS
This file aims to provide more details concerning the QTL x Year interaction of the resistance of Durum Wheat to WSSMV. We are going to study SS. We will use QTLRel
library(QTLRel)## R/QTLRel is loaded
library(xtable)Introduction
We are going to study the effects of markers trough years and pop progressively.
Load data
I load the genotyping matrix first.
genotype<-read.table("/Users/holtz/Dropbox/Publi_Mosaique/DATA/DATA/GROUPED/genotypage.csv", sep = ";" , header = F, na.strings = "-")
genotype=as.matrix(genotype)
colnames(genotype)=genotype[1,]
genotype=as.data.frame(genotype[-1 , ])
names(genotype)[1]<-"geno"
print("--- Your genotyping matrix looks correct. Dimension of the matrix are :")## [1] "--- Your genotyping matrix looks correct. Dimension of the matrix are :"
print(dim(genotype))## [1] 348 7342
# I copy this matrix 2 times, since I read 2012 and 2015 together.
rownames(genotype)=genotype[,1]
genotype=genotype[,-1]
a=genotype ; rownames(a)=paste(rownames(a),"2012",sep="_")
b=genotype ; rownames(b)=paste(rownames(b),"2015",sep="_")
genotype=rbind(a,b)Then the genetic map:
map <- read.table("/Users/holtz/Dropbox/Publi_Mosaique/DATA/DATA/genetic_map.txt" , header=T , dec = ".", na.strings = "-" , check.names=F)
colnames(map) <- c("LG", "marqueur", "Distance","group_physique","Posi_physique")
rownames(map) <- map$marqueur
map$LG <- as.factor(map$LG)
print("--- Your genetic map looks correct. Dimension of the map are :")## [1] "--- Your genetic map looks correct. Dimension of the map are :"
print(dim(map))## [1] 8568 5
map=map[ , c(2,1,3,5)]
colnames(map)=c("snp","chr", "dist", "phyPos")And finally the phenotyping matrix
BLUP<-read.table("/Users/holtz/Dropbox/Publi_Mosaique/DATA/DATA/GROUPED/phenotypage.csv", header = TRUE, sep=";")
colnames(BLUP)[1]="geno"
print("--- Your Phenotyping matrix looks correct. Dimension of the matrix are :")## [1] "--- Your Phenotyping matrix looks correct. Dimension of the matrix are :"
print(dim(BLUP))## [1] 396 6
# Fichier de phénotypage modifié, il va falloir mettre la SS de 2012 et 2015 ensemble, avec une colonne année.
a=BLUP[, c(1,2)] ; a$year="2012" ; colnames(a)=c("geno", "SS_blup_AR1","year")
#a[,2]=a[,2]/sqrt( mean(c(0.66,0.83)) )
#a[,2]=ifelse(substr(a$geno, 1,2)=="TT", a[,2]/sqrt(0.66) , a[,2]/sqrt(0.83))
b=BLUP[, c(1,4)] ; b$year="2015" ; colnames(b)=c("geno", "SS_blup_AR1", "year")
#b[,2]=b[,2]/sqrt( mean(c(1.2,1.21)) )
#b[,2]=ifelse(substr(b$geno, 1,2)=="TT", b[,2]/sqrt(1.2) , a[,2]/sqrt(1.21))
BLUP=rbind(a,b)
rownames(BLUP)=paste( BLUP[,1], BLUP$year,sep="_")
BLUP=BLUP[,-1]
# Note: On peut garder les blups tels quels / ou les pondéré par la variance génet de chaque année moyennée sur les 2 pops / ou par la variance génet de chaque année et chaque pop.Prepare data
We need to have genotype and phenotype in the same order.
And I add a “pop” column in the phenotyping matrix:
Y=na.omit(BLUP)
Y=Y[which(rownames(Y)%in%rownames(genotype)) , ]
Y$pop=substr(rownames(Y),1,2)
genotype=genotype[which(rownames(genotype)%in%rownames(Y)) , ]
genotype=genotype[ match(rownames(Y),rownames(genotype)) , ]Impute missing data
# missing data
set.seed(123)
XNNA=genotype
my_fun=function(x){length(x[x=="A" & !is.na(x) ])/length(!is.na(x)) }
prop=apply(XNNA , 2 , my_fun)
for(i in c(1:ncol(XNNA))){
aa=XNNA[,i][is.na(XNNA[,i])]
bb=rbinom(length(aa),1,prob=prop[i])
XNNA[,i][is.na(XNNA[,i])]=c("B","A")[bb+1]
}
# Change "A" to "AA"
XNNA=as.matrix(XNNA)
XNNA[which(XNNA=="A")]<-"AA"
XNNA[which(XNNA=="B")]<-"BB"test1: DS - 2012
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which(Y$pop=="TT" & Y$year=="2012") , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which (substr(rownames(XNNA),1,2)=="TT" & grepl("2012" , rownames(XNNA))) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 161 3544
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DS_2012=AAtest2: DS - 2015
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which(Y$pop=="TT" & Y$year=="2015") , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which (substr(rownames(XNNA),1,2)=="TT" & grepl("2015" , rownames(XNNA))) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 161 3544
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DS_2015=AAtest3: DL - 2012
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which(Y$pop=="BX" & Y$year=="2012") , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which (substr(rownames(XNNA),1,2)=="BX" & grepl("2012" , rownames(XNNA))) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 181 5851
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DL_2012=AAtest4: DL - 2015
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which(Y$pop=="BX" & Y$year=="2015") , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which (substr(rownames(XNNA),1,2)=="BX" & grepl("2015" , rownames(XNNA))) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 185 5851
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DL_2015=AAtest5: DS AND DL - 2012
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which( Y$year=="2012") , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which ( grepl("2012" , rownames(XNNA))) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 342 7341
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, x=Y_tmp$pop, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$pop, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DSDL_2012=AAtest6: DS AND DL - 2015
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which( Y$year=="2015") , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which ( grepl("2015" , rownames(XNNA))) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 346 7341
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, x=Y_tmp$pop, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$pop, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DSDL_2012=AAtest7: DS - 2012 & 2015
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which(Y$pop=="TT" ) , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which (substr(rownames(XNNA),1,2)=="TT" ) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 322 3544
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DS_2012_2015=AAtest8: DL - 2012 & 2015
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y[which(Y$pop=="BX" ) , ]
# Select the corresponding genotyping data?
XNNA_tmp=XNNA[ which (substr(rownames(XNNA),1,2)=="BX" ) , ]
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 366 5851
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DL_2012_2015=AAtest9: DS AND DL - 2012 & 2015
QTL detection with QTL rel for DS only in 2012 only
# Select the corresponding phenotype data
Y_tmp=Y
# Select the corresponding genotyping data?
XNNA_tmp=XNNA
XNNA_tmp=XNNA_tmp [ , which(apply(XNNA_tmp , 2 , function(x){length(unique(x))} )==2) ]
dim(XNNA_tmp)## [1] 688 7341
# Kinship matrix
K<-genMatrix(XNNA_tmp)
# I = identity matrix
I<-diag(nrow(Y_tmp))
# Variance components
mod1<-estVC(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, v=list(AA=K$AA,DD=NULL,HH=NULL,AD=NULL,MH=NULL,EE=I))
# Test marker per marker
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, vc=mod1, gdat=XNNA_tmp ,test="Chisq")Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(llk.hk1$p) , llk.hk1$p), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$llk.hk1.p) , ]
# And plot it
plot(-log10(AA$llk.hk1.p) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DL_2012_2015=AAMore info on the 2 weird significant QTLs
bilan<-data.frame(marqueurs=names(llk.hk1$p),LOD=-log10(llk.hk1$p), r2=llk.hk1$v)
bilan=merge(map,bilan, by.x=1 , by.y=1, all.y=T)
bilan=bilan[order(bilan$chr, bilan$dist) , ]
a=bilan[which(bilan$chr=="2B") , ]
a=a[which(a$LOD==max(a$LOD, na.rm=T)) , ]
b=bilan[which(bilan$chr=="5B") , ]
b=b[which(b$LOD==max(b$LOD, na.rm=T)) , ]
bilan=rbind(a,b)
head(bilan)## snp chr dist phyPos LOD
## 377 Cluster_10786|Contig2|original@200 2B 117.8 44285664 4.292400
## 6270 Cluster_8710|Contig2|original@224 5B 99.8 183184647 4.806105
## r2
## 377 2.338593
## 6270 2.650581
Distortion de ségrégation?
tmp=data.frame(all=genotype[ , "Cluster_10786|Contig2|original@200"], pop=substr(rownames(genotype), 1 , 2))
table( tmp )## pop
## all BX TT
## A 366 136
## B 0 182
tmp=data.frame(all=genotype[ , "Cluster_8710|Contig2|original@224"], pop=substr(rownames(genotype), 1 , 2))
table( tmp )## pop
## all BX TT
## A 133 130
## B 219 186
test10: DS AND DL - 2012 & 2015 - interaction effect
QTL detection with QTL rel for DS only in 2012 only
# Beginning is the same as test 9
# Test marker per marker with and without interaction
llk.hk1 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, vc=mod1, gdat=XNNA_tmp ,test="None")
llk.hk2 <- scanOne(y=Y_tmp$SS_blup_AR1, x=Y_tmp$year, intcovar=Y$year, vc=mod1, gdat=XNNA_tmp ,test="None")
# Find pvalue of each marker for interaction significance
diff_deviance <- llk.hk2$p - llk.hk1$p
pvaluesexpected=1-pchisq(diff_deviance, 1)Let’s observe results.
# Merge LODs with the genetic map
AA=merge(map,data.frame( names(pvaluesexpected) , pvaluesexpected), by.x=1 , by.y=1, all.x=T)
AA=AA[order(AA$chr, AA$dist) , ]
AA=AA[!is.na(AA$pvaluesexpected) , ]
# And plot it
plot(-log10(AA$pvaluesexpected) , pch=20 , col=as.numeric(AA$chr) , cex=1.3, xaxt="n" , ylim=c(0,4) )
abline(h=-log10(0.05), col="grey", lwd=1.5)
abline(h=3.61, col="grey", lwd=1.5)
num=seq(1,nrow(AA))
num=aggregate(num, by=list(AA$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])# Save as the object
res_DSDL_2012_2015_inter=AAShow the 3 interaction effects (marker details)
a=res_DSDL_2012_2015_inter[which(res_DSDL_2012_2015_inter$chr=="2A") , ]
a=a[a$pvaluesexpected==min(a$pvaluesexpected) , ]
b=res_DSDL_2012_2015_inter[which(res_DSDL_2012_2015_inter$chr=="7A") , ]
b=b[b$pvaluesexpected==min(b$pvaluesexpected) , ]
c=res_DSDL_2012_2015_inter[which(res_DSDL_2012_2015_inter$chr=="7B") , ]
c=c[c$pvaluesexpected==min(c$pvaluesexpected) , ]
signif_inter=rbind(a,b,c)| snp | chr | dist | phyPos | pvaluesexpected | |
|---|---|---|---|---|---|
| 1526 | Cluster_14091|Contig2|original@276 | 2A | 218.80 | 240000045 | 0.00 |
| 2143 | Cluster_16523|Contig1|complementarySeq@373 | 7A | 119.70 | 43724948 | 0.01 |
| 2161 | Cluster_1658|Contig1|likelySeq@900 | 7B | 43.80 | NA | 0.00 |
LOD Threshold
#===== FORMAT R-QTL
# fichier genfile
donnees <- merge(genotype, Y , by="row.names")
donnees=donnees[ , -which(colnames(donnees)%in%c("year","pop")) ]
indiv <- donnees[,1]
tmarq <- t(donnees[, -ncol(donnees)])[-1, ]
colnames(tmarq) <- donnees[,1]
tmarq <- data.frame(marqueur = rownames(tmarq), tmarq)
fich <- merge(map, tmarq, by.x=1 , by.y=1 ,sort = FALSE)
tfich <- t(fich)
tfich <- data.frame(ID = c("ID", "", "", rownames(tfich)[4:length(rownames(tfich))]), tfich)
write.table(tfich, file = "tmpgen.csv", row.names = FALSE, sep = ";", col.names = FALSE,quote = FALSE)
# fichier phefile
phen <- donnees[, c(1, ncol(donnees) )]
colnames(phen)[1] <- "ID"
write.table(phen, file = "tmpphe.csv", row.names = FALSE, sep = ";", col.names = TRUE, quote = FALSE)
# ===== CALCUL DU THRESHOLD VIA rQTL
library(qtl)
set.seed(123)
donnees <- read.cross("csvs", genfile = "tmpgen.csv", phefile = "tmpphe.csv" , sep=";")
donnees <- convert2riself(donnees)
data=calc.genoprob(donnees)
res_permut <- scanone(cross = data, pheno.col = 2, method = "hk", n.perm = 1000)
a=summary( res_permut, alpha = 0.05 )
print(my_carac)
print(a)–> 3.62. Did not change.